#include "triangle.h"

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <limits>
// Global variable... Bad bad Programmer, no donut for you!

boost::mt19937 mersenne_twister;
boost::uniform_real<> real_dist(0.0,1.0);
boost::variate_generator<boost::mt19937&, boost::uniform_real<> > random_number(mersenne_twister, real_dist);

namespace t2c2 {
	
Triangle::Triangle(Material* _material,const Vector3D& _vertex1,const Vector3D& _vertex2,const Vector3D& _vertex3) : Patch(_material) {
	vertices[0] = _vertex1;
	vertices[1] = _vertex2;
	vertices[2] = _vertex3;
}

Triangle::Triangle(const Vector3D& _vertex1,const Vector3D& _vertex2,const Vector3D& _vertex3) : Patch() {
	vertices[0] = _vertex1;
	vertices[1] = _vertex2;
	vertices[2] = _vertex3;
}

Triangle::Triangle(void) : Patch() {
	
}

Triangle::~Triangle(void){

}

Vector3D Triangle::get_vertex(unsigned short _index) const {
	return vertices[_index];
}

void Triangle::set_vertex(unsigned short _index,const Vector3D& _vertex) {
	vertices[_index] = _vertex;
}

double Triangle::area(void) const {
	const Vector3D edge1 = vertices[1] - vertices[0];
	const Vector3D edge2 = vertices[2] - vertices[0];
	return 0.5*(edge1.cross_product(edge2)).length();
}

Vector3D Triangle::get_normal(const Vector3D&) const {
	const Vector3D edge1 = vertices[1] - vertices[0];
	const Vector3D edge2 = vertices[2] - vertices[0];
	return (edge1.cross_product(edge2)).normalize();
}

Vector3D Triangle::get_sample_point(void) const {
	double u;
	double v;
	do {
		u = random_number();
		v = random_number();
	} while((u + v) > 1.0);

	//Barycentric coordinates
	Vector3D sample_point = vertices[0]*(1.0 - u - v) + vertices[1]*u + vertices[2]*v;
	
	return sample_point;
}
double Triangle::intersection(const Ray& _ray) {
  // Tomas Möller and Ben Trumbore ray-triangle intersection

	std::numeric_limits<double> double_limit;
	const double infinity = double_limit.infinity();
	const double epsilon = double_limit.epsilon();
  //const double epsilon = 0.000001;

	Vector3D intersection_point;

	Vector3D edge1 = vertices[1] - vertices[0];
	Vector3D edge2 = vertices[2] - vertices[0];

	Vector3D pvec = _ray.get_direction().cross_product(edge2);

	double det = edge1.dot_product(pvec);
	if(get_normal(Vector3D()).dot_product(_ray.get_direction()) <  0.0) return -1;
  // Culling test
	if(det > -epsilon && det < epsilon) return -1;
	double inv_det = 1.0 / det;

	Vector3D tvec = _ray.get_origin() - vertices[0];

	double u = tvec.dot_product(pvec) * inv_det;

	if(u < 0.0 || u > 1.0) return -1;

	Vector3D qvec = tvec.cross_product(edge1);

	double v = _ray.get_direction().dot_product(qvec) * inv_det;

	if(v < 0.0 || u + v > 1.0) return -1;

	double t = edge2.dot_product(qvec) * inv_det;
	if(t < epsilon) return -1;

	intersection_point = vertices[0]*(1.0 - u - v) + vertices[1]*u + vertices[2]*v;
	return t;                                                                        
}

}
